#PSRM Replication Files
#
#The Population Ecology of Interest Groups and Counter-Mobilization: 
#Reproductive Rights Organizations in the United States, 1920-1985.
#
#Tristan M. Hightower (2023)

#Change working directory for data loading

setwd("~/Desktop/PSRM Replication")

#Load libraries
library(readr)
library(imputeTS)
library(readxl)
library(MASS)
library(stargazer)
library(tidyverse)
library(growthcurver)

#Load Data
groups <- read_csv("PSRM_Replication.csv")

#Descriptive Figure 1

plot(groups$Year,groups$Choice_density, type = "l", xlab = "Year",
     ylab = "Density", xlim = c(1930,1985), ylim = c(0,70))

lines(groups$Year, groups$Life_density, lty = 3)

legend("topleft",
       legend=c("Pro-Choice", "Pro-Life"),
       lty=c(1,3))


abline(v=1973, lty = 5)
text(1967.5,60,"Roe v. Wade", cex = .8)


#Run Models
m1<-glm.nb(groups$Life_founding~groups$Life_density+groups$Life_densitysq+
             groups$Choice_density+groups$Choice_densitysq+groups$Abortion_ratio+
             groups$M_BA+groups$W_BA)

m2<-glm.nb(groups$Choice_founding~groups$Choice_density+groups$Choice_densitysq+
             groups$Life_density+groups$Life_densitysq+groups$Abortion_ratio+groups$M_BA+groups$W_BA)


#Create Table 1
  
stargazer(m1,m2, type = "text",
          column.labels = c("Anti-Abortion","Abortion-Rights"),
          dep.var.labels.include = F, dep.var.caption = "DV: Group Formation",
          covariate.labels = c("Anti-Abortion", "Anti-Abortion Sq.",
                               "Abortion Rights", "Abortion Rights Sq.",
                               "Abortion Ratio", "Male BA Degrees","Female BA Degrees"),
          notes = "Negative binomial regression with standard errors in parentheses.")



#Predictions

predictions <- data.frame("Choice_density" = 0:65, "Life_density" = 0:65, "Abortion_ratio"= mean(groups$Abortion_ratio))
predictions$Life_densitysq<-(predictions$Life_density)^2
predictions['Choice_foundings']<-NA
predictions$Choice_densitysq<-(predictions$Choice_density)^2

predictions['Life_foundings']<-NA

predictions$Choice_foundings<-predict(m2,newdata=predictions)
predictions$Life_foundings<-predict(m1,newdata=predictions)

#Plotting Predictions


#Figure 2, publication
ggplot(predictions, aes(Life_density,Choice_foundings)) +
  geom_smooth(level = .9, method = "lm") +
  xlab("Density of Anti-Abortion Groups")+
  ylab("Predicted Abortion Rights Formation")+
  labs(title = "Pred. Abortion-Rights Foundings Relative to Anti-Abortion Density") +
  coord_cartesian(xlim = c(0,50))

#Figure 3, publication
ggplot(predictions, aes(Life_density, Life_foundings)) +
  geom_smooth(level = .9, method = "lm") +
  xlab("Density of Abortion-Rights Groups")+
  ylab("Predicted Anti-Abortion Foundings")+
  labs(title = "Pred. Anti-Abortion Foundings Relative to Abortion-Rights Density") 


###Appendices###

#Creating figures 4 and 5 using growthcurver 'Summarize_Growth' function

#creating blank dataset into 'd2'
groups$time<-c(1:66)
d2<-groups$time
d2<-as.data.frame(d2)
colnames(d2)<-"time"
d2$choice<-groups$Choice_density
d2$life<-groups$Life_density

#running model for abortion rights
gc_fit_choice<-SummarizeGrowth(d2$time,d2$choice)

#Creating Figure 4
plot(gc_fit_choice, xlab = "Time (Years)", ylab = "Density", main = "Abortion-Rights Growth Curve")

#running model for anti-abortion
gc_fit_life<-SummarizeGrowth(d2$time,d2$life)


#Creating Figure 5
plot(gc_fit_life, ylab="Density", xlab = "Time (Years)", main = "Anti-Abortion Growth Curve")



#Creating Table 2
groups2<-groups[groups$Year>1966,]
m3<-glm.nb(groups2$Life_founding~groups2$Life_density+groups2$Life_densitysq+groups2$Abortion_ratio)
m4<-glm.nb(groups2$Life_founding~groups2$Choice_density+groups2$Choice_densitysq+groups2$Abortion_ratio)

stargazer(m3,m4,dep.var.caption = "DV: Anti-Abortion Group Formation",dep.var.labels.include = F,
          type = "text", covariate.labels = c("Anti-Abortion", "Anti-Abortion Sq.", "Abortion Rights", 
                                              "Abortion Rights Sq.", "Abortion Ratio"),
          omit = c("Constant"), omit.stat = "theta", notes = "Negative binomial regression with SEs in parentheses.")
